Simulating dark-field X-ray microscopy images with wavefront propagation techniques

The simulation of a dark-field X-ray microscopy experiment using wavefront propagation techniques and numerical integration of the Takagi–Taupin equations is shown. The approach is validated by comparing with measurements of a near-perfect diamond crystal containing a single stacking-fault defect.


Introduction
Dark-field X-ray microscopy (Simons et al., 2015) (DFXM) is a full-field X-ray imaging technique similar to X-ray topography (Berg, 1931;Lang, 1957). However, unlike the latter, DFXM utilizes an X-ray objective lens placed between the sample and the detector to create a magnified image and can therefore achieve a spatial resolution better than the detector pixel size. The spatial bandwidth is limited by the small numerical aperture (NA ' 10 À3 ) of this objective lens. Compared with classical X-ray topography, this provides angular resolution of the scattered beam direction (Poulsen et al., 2017), which makes it possible to quantitatively measure strains and rotations of the crystal lattice by sequentially collecting images while rotating the sample and moving the objective lens.
Traditionally, the quantitative analysis of DFXM data and the theoretical description (Poulsen et al., 2017 of the method have relied on strong approximations. Most important is the kinematical approximation, which is to omit multiple scattering effects, and which holds for small and for highly deformed crystals. Another approximation, which is built into the geometric optics treatment of Bragg scattering, is that infinitesimal sub-volumes in the sample scatter according to the Bragg law for a perfect infinite crystal, and that the intensities scattered from such sub-volumes add together incoherently. In some cases (e.g. for near-perfect crystallites and small defect structures), however, such approximations are not valid.
Here we present a method for simulating DFXM images based on wavefront propagation combined with a framework that treats multiple scattering events (known as the dynamical theory of X-ray diffraction) (Authier, 2001). This method is able to handle the effects of coherence, dynamical scattering, aberrations of the objective lens and detector imperfections, and as such it is a more realistic model of DFXM for nearperfect crystals than the geometric optics model. This will be useful for understanding the type of contrast observed in DFXM images and can aid in experimental planning and data analysis.
Furthermore, a number of advanced approaches to DFXM have been suggested in the literature and tested by various authors [e.g. magnified topo-tomography , confocal Bragg microscopy and tele-ptychography (Pedersen et al., 2020), Fourier ptychographic DFXM (Carlsen et al., 2022)]. To date, the theoretical models used in these examples rely on idealized models of the underlying physics. Again, a more realistic forward model would be useful for testing the limits of their applicability.
In this paper we discuss the steps that are required for such a simulation and present an implementation and some initial results. We compare simulated and experimental data for a near-perfect single-crystal diamond containing only one single stacking fault in the imaged field of view (FOV).

Method
Recent developments in DFXM are moving towards studying the dynamics of isolated defects, such as dislocations, domain walls and acoustic waves in near-perfect single crystals Holstad et al., 2021). Here, dynamical effects cannot be ignored. In these cases, the analysis of DFXM data has often relied on the weak beam approximation, where it is assumed that even highly perfect crystals scatter approximately kinematically when the sample is rotated to the tails of the rocking curve. At this position, the bulk of the crystal does not scatter strongly. Instead, only small strained volumes near defects and surfaces contribute to the scattered intensity.
Much has been published on the subject of simulating X-ray topography images (Taupin, 1964;Authier, 1968;Epelboin & Soyer, 1985). These methods are applicable in our case, but a few extra precautions must be taken due to our need for high quantitative accuracy, and the added complication caused by the objective lens. There exist useful solutions (Chubar et al., 2013;Pedersen et al., 2018;Celestre et al., 2020) for simulating synchrotron sources and X-ray optical components of the kinds applied in a DFXM instrument with coherent wavefront methods. A full simulation based on coherent wavefronts is therefore within reach by combining these established methods.
In dark-field transmission electron microscopy (DFTEM), simulating images with dynamical diffraction effects is done routinely. Although superficially similar to DFXM, the methods cannot be transferred from one technique to the other. The most common methods in DFTEM are either Bloch-wave methods that rely on the column approximation which cannot be applied to DFXM where the scattering angles often exceed 10 , or multi-slice methods that require the atomic structure of the sample to be well sampled, which is not feasible for DFXM where sample sizes often exceed 100 mm (Kirkland, 1998).
X-ray topography methods, on the other hand, rely on the two-beam approximation (Taupin, 1964), which cannot be applied to electron microscopy where potentially hundreds of reflections contribute significantly to the images.
Here we apply the formalism of the Takagi-Taupin equations, which make use of the two-beam approximation to simplify the scattering problem. This allows us to treat Bragg scattering in an averaged way, removing the requirement to over-sample the unit cell. Instead, the sampling is limited by the divergence of the incident X-rays and the size of features in the sample. Propagation from the sample to the detector is handled by established paraxial Fourier optics methods.
We split the simulation flow into five steps as follows: (i) Beam: calculate the amplitudes of the modulated waves of the beam incident on the sample.
(ii) Sample: calculate the complex scattering function throughout the sample crystal.
(iii) Integration: numerically integrate the Takagi-Taupin equations to get the complex amplitudes of the scattered beam on the exit surface of the sample.
(iv) Propagation: propagate the scattered complex amplitudes through the imaging optics to the detector.
(v) Detector: interpolate the propagated field to the detector pixels and account for detector characteristics.
A schematic drawing of the steps and flow of the simulation is shown in Fig. 1.

Defining the beam
We work in a paraxial wave-optics formalism and use a coherent mode decomposition to describe the state of the X-ray beam at a given plane. This is given by a number of modes, labeled with p: E (p) (x, y; z), where E is the amplitude function of a modulated plane wave:  E ¼ <½EðrÞ expðik 0 Á r À i!tÞp p, where r = (x, y, z) is the spatial position, k 0 is the modes' wavevector, h -! is the photon energy andp p is the polarization vector.
If a mode is known on a plane z = 0, then the mode's amplitude on another plane can be found by applying a (linear) coherent wavefront propagator: E ðpÞ ðx; y; zÞ ¼ P 0!z ½E ðpÞ ðx; y; 0Þ: We assume that the radiation is beam-like, i.e. that it has a limited extent in two dimensions in real space and in all dimensions in reciprocal space. The reciprocal-space distribution is around a central wavevector k 0 , which defines the nominal direction and wavelength of the incident beam.

Defining the sample
In the Takagi-Taupin (Takagi, 1962(Takagi, , 1969Taupin, 1964) approach to X-ray scattering in the two-beam approximation, the only quantities of interest are the average electric susceptibility of the crystal, 0 , and the two scattering constants h and h that describe the cross section for scattering and back-scattering, respectively. These may be spatially modified by a displacement field u(r) of the strained crystal: Here Q is the reciprocal-lattice vector of the given reflection in an undeformed reference lattice. This dependence on the displacement field explains the high sensitivity to small strains. h and h can often be considered constants, but in samples with twin boundaries there may be discontinuous jumps in the values of this function, which explains contrast observed at presumably strain-free inversion twin domain boundaries in polar materials (Klapper, 1987).
In this paper we focus on slab-shaped single crystals, i.e. crystals with two parallel surfaces and infinite size in the orthogonal directions (see Fig. 2). The normal of the exit surface is calledn n.
We utilize a discrete representation of the sample structure on an orthogonal grid defined by the three directionsx x,ŷ y and z z ¼n n, and corresponding step sizes d x , d y and d z . It is necessary that the surface normal of the crystal slab is parallel to thê z z axis, but we make no requirement on the last free rotation. The number of grid points in each direction will be labeled N x , N y and N z , giving a total size of the simulated domain of L x = d x (N x À 1), L y = d y (N y À 1), L = d z (N z À 1). L is the thickness of the simulated crystal slab.
The complex value of these scattering functions needs to be known with high resolution. For simple test cases, where the displacement field is given by an analytical expression, this is not a point of concern. If, however, the displacement field is generated by a numerical simulation, then it needs to be computed with sufficient resolution to at least match the resolution of the experiment ('50 nm) throughout the volume of the sample. The size of the sample can be up to several hundred micrometres.
For highly deformed samples, the scattering function contains a phase factor of the shape exp(iQ Á u). In order to limit the phase variation between adjacent voxels to less than 2, step sizes must be below |ru|/(2|Q|), which means small steps must be used for highly deformed samples. In samples with nanometre-sized domains or other small structures, all structural features must be resolved.

Integrating the Takagi-Taupin equations
Scattering of X-rays by a deformed crystal is treated in the formalism of the Takagi-Taupin equations (TTEs). In formulating the TTEs, there is an arbitrary choice of the vector k 0 , which leads to slightly different versions of the final equations. The ones chosen here are referred to as the symmetrical TTEs (Vartanyants & Robinson, 2001) that arise when k 0 is taken to be the vacuum wavevector of the incident beam as opposed to the more conventional choice of the refracted wavevector used in the original publications: where ¼ 2 sinð2Þ is a measure of the misorientation away from the vacuum Bragg condition, is the rocking angle and C is the polarization factor. E ðpÞ 0 are the modes of the plane-wave decomposition of the incident beam and E ðpÞ p are the corresponding modes of the scattered beam, which is given relative to the wavevector k h = k 0 + Q.
The TTEs are typically solved using a finite-difference integration scheme. A number of different algorithms have been used in the literature with the pioneering work done by Taupin (1964), Authier (1968), Epelboin & Soyer (1985). In this paper we use a novel, more flexible scheme which is documented by . In contrast to other established methods, this method is able to operate on an orthogonal grid representation of the sample structure, which is achieved by implicitly utilizing Fourier interpolation. This, however, causes a less efficient use of computer resources.
The geometry of the problem is set by the shape of the incident beam, the vectors k 0 and Q, as well as the choice of a computational grid. The plane spanned by the two vectors k 0 and Q is called the scattering plane.
When the scattering plane is normal to theŷ y axis, the TTEs decouple into a set of 2D problems that can be solved slice-byslice. If we are free to choose the orientation of the computational grid, we can always choose a geometry where this becomes the case.
When Q is parallel to the surface of the crystal, we say that we have a symmetric Laue geometry. In this case we can choose an orientation of the computational grid where Qjjx x, where the TTEs take a particularly simple form.
In order to solve the TTEs, we impose zero Dirichlet boundary conditions in the two transverse dimensions, x and y. These boundary conditions require the sample grid to be large enough to fit the entire Borrmann triangle extending from every point where the incident beam amplitude is non-zero. This is fulfilled if the non-zero part of the amplitude is fully contained in the rectangle defined by (see Fig. 2 where L x , L y are the lengths of the simulated domain in the x and y directions, and L is the thickness of the crystal. For a given incident beam and Q, this sets a minimum on the size of the simulated domain in the two transverse directions. With the sample structure and the boundary conditions given, the TTEs constitute an initial value problem, where the initial value is the amplitude of the incident X-rays on the z = 0 surface. This can be solved by an appropriate finite-difference scheme to yield the amplitudes of both the transmitted and scattered beams on the exit surface z = L.

Propagating through the optics
In DFXM, the scattered beam at the exit surface of the sample is imaged onto a detector using an objective lens in a magnifying geometry. The propagation through the lens is challenging to simulate due to the inherent thick-lens behavior of the compound refractive lens (CRL) which is typically used as the objective lens. For this, we use a computational approach where each lens in the CRL is treated as a thin lens and a paraxial FFT (fast Fourier transform) propagator is used to propagate the wavefront between each lens.
To do this, we make use of a method for propagating wavefronts with rapidly varying quadratic phases that come from the transmission function of the lenses that treats the quadratic component analytically (Ozaktas et al., 1996;Chubar & Celestre, 2019). To simulate an aberrated CRL, a separate aberration function is included at each lenslet in the CRL.
These methods have previously been used to model thick CRLs like the one used in this study (Pedersen et al., 2018;Celestre et al., 2020).
It is useful to introduce a new optical axis (ẑ z i ) aligned with the average wavevector of the scattered beam, k h . To bridge the gap between these two coordinate systems, we project the values from the exit surface of the crystal to the z = 0 plane of the new coordinate system.
The inherent near-field nature of the imaging geometry, which is determined by the fact that the FOV is as large as the aperture of the objective lens, is handled by first multiplying the field with the near-field phase factors: The need to over-sample this function can effectively limit the size of the FOV that is possible to simulate for a given pixel size.

Detector characteristics
With the mode amplitudes on the detector plane given, we now need to interpolate these values to the detector pixels and incoherently sum over the modes of the incident beam.
If the purpose of the simulation is to estimate the resolution or to create a data set to be used to test the data analysis procedures, one should remember to include non-ideal behavior introduced at the detector. The most important effects are the incoherent point spread of the detector, background signal, counting noise and non-linear response.
The detector used here is an indirect detector composed of a scintillator crystal (25 mm-thick gadolinium gallium garnet), optical microscope (Mitutoyu M Plan Apo 10Â, NA = 0.28 and tube lens) and pco.edge 5.5 sCMOS camera with pixel size 6.5 mm. Contributions to the incoherent point spread function arise from the diffraction limit due to the finite NA of the optical microscope, the finite thickness of the scintillator (especially when it is larger than the depth-of-focus of the optical microscope), scattering within the scintillator and the optical microscope, and aberrations.
Counting statistics/readout noise can be the critical factor when imaging small grains or weak reflections. The non-linear response (saturation) might be quite important for perfect crystals, as the images have interesting features over a very large dynamic range.

Comparison with experiment
To test our approach, we simulate a section-topography-type experiment with a near-perfect single-crystal diamond slab of thickness 300 mm containing a single stacking-fault defect. The sample is imaged in a symmetric Laue geometry in a {111} reflection with [110] entrance and exit surfaces.
The investigated defect is a stacking fault, which arises by the addition or removal of a single close-packed plane of atoms in the face-centered cubic (f.c.c.) parent lattice of the diamond. The fault vector is of the family b s:f: ¼ ð1=3Þh111i which is not a translational symmetry of the f.c.c. lattice. These planar defects are bounded by the surfaces of the crystal and Frank-type partial dislocations (Frank, 1951;Kowalski et al., 1989). In the described experiment, the edges of this defect lie outside the Borrmann triangle and the effects of the strain originating at the edges can be ignored. This allows us to treat the stacking fault as an infinite planar defect [on one side u(r) = 0, on the other u(r) = b s.f. ]. In the Takagi-Taupin description, the stacking fault thus becomes a discrete jump in the phase of the scattering function of magnitude 2[hk'] Â b s.f. = AE2/3 when imaged using a [111] reflection not orthogonal to the stacking-fault normal (Klapper, 1987). A sketch of the sample geometry is shown in Fig. 3.
This constitutes a good test sample, as the defect (aside from the 12 possible orientations of the defect) only has a single continuous degree of freedom: the position of the stacking fault along the normal. Furthermore, diamond is one of a few materials where macroscopic crystals with very low defect density are available. Provided our method captures all the relevant physics, we should therefore be able to perfectly recreate the experimental data.
The dynamical scattering patterns produced by isolated stacking faults in diamonds have previously been studied in detail by classical X-ray topography (Kowalski et al., 1989).
Experiments were carried out at the ESRF dark-field X-ray microscopy beamline, ID06-HXM (Kutsal et al., 2019), after the EBS (Extremely Brilliant Source) upgrade. A Si (111) Bragg-Bragg double-crystal monochromator selected X-rays with photon energy 17 keV from the undulator source. The spot size of the beam on the condenser lens is limited by a slit in the vertical direction to 0.2 mm. Based on parameters of the X-ray source given by the ESRF, we estimate a vertical coherence length of 340 mm at the distance of 50 m from the source to the condenser lens, which is significantly larger than the aperture of the condenser lens. The horizontal coherence length is estimated to be 40 mm, which is smaller than the FOV but much larger than the extent of the coherent point spread function of the objective lens. The objective lens consists of 70 individual biconcave Be lenses of apex curvature 50 mm. The CRL is further modified with a 100 mm square aperture after the last lens To describe the incident X-rays, we use only one coherent mode for each wavelength and 41 different wavelengths covering a relative energy spread of 6 Â 10 À4 in total. For each energy component, the amplitude at the sample position is the Fourier transform of the complex transmission function of an aberration-free 1D condenser lens with a Gaussian-shaped absorbing aperture and a hard cut-off at 140 mrad. Consequently, the angular spectrum in the transverse direction is a top-hat profile with a small Gaussian smoothing at the edges, whereas the spatial profile is a diffraction-limited focal spot with some ringing due to the hard cut-off (see Fig. 4). We omit the use of several transverse mutually incoherent modes. The effects of including these would be small, since the effective source size is small compared with the NA of the objective lens.
In the experimental realization we make the observations in a laboratory frame defined byx x lab ,ŷ y lab andẑ z lab , whereẑ z is parallel to k 0 . In this experiment the scattering is vertical in the laboratory frame, meaning thatŷ y lab is normal to the scattering plane. We choose the integration grid such thatŷ y ¼ŷ y lab .
The simulation used a grid of 2048 Â 2048 Â 3001 points with step sizes of 60, 60 and 100 nm, respectively. This matches the 300 mm thickness of the sample and the '100 mm FOV in the y direction. The large size in the x direction was necessary to satisfy the constraints of equation (5). Step sizes in x and y directions are chosen to be larger than the resolution of the final images, and the step size in z is chosen such that the integration error of the finite-difference scheme is lower than the noise level of the final images. The execution time is dominated by the integration of the TTEs, which took 3.5 h per energy point running on a single core. The simulation was parallelized over the modes.
The computation time was dominated by the integration of the TTEs, which involved '10 000 1D Fourier transforms of the 2048 Â 2048 arrays used in this example compared with just '150 2D Fourier transforms for the simulation of the optics. The use of computer resources scales linearly with the number of modes, so if a transversely incoherent model was to be used, the needed resources would increase by a factor of the number of modes. In the study presented here computa-   Vertical profile in (a) reciprocal and (b) real space of the beam used in the simulation of this paper. tion time was not an issue, but if high-throughput simulation was needed, large performance gains could likely be achieved by applying a more efficient algorithm for the integration of the TTEs (Epelboin & Soyer, 1985). The integration of the TTEs can also be parallelized into separate 2D problems for each vertical slice of the sample and should easily scale with available computer resources. If the thick-lens effects are not critically important, the CRL may well be approximated by a single convolution step (Pedersen et al., 2018). Furthermore, the large number of energy modes used in this example was only necessary to ensure proper sampling of the highfrequency features used as lens aberration and avoid aliasing problems in Fig. 9. If these are not important, a much smaller number of energy points would be sufficient.

Near-field measurements
The microscope used for the experiment provides the possibility to additionally carry out traditional X-ray topography by placing a detector closely behind the sample (40 mm in the examples shown here) without using the objective lens. This is useful for alignment and for lowresolution characterization of the sample.
For a single coherent mode [ Fig. 5(b)] this propagation distance causes recognizable Fresnel-diffraction fringes around sharp features in the scattered field. In a polychromatic simulation [Fig. 5(c)], these fringes are blurred out. The difference in wavelength (which causes a small difference in the free-space propagator) is not sufficient to explain this blurring. Rather, the broadening is caused by the slight difference in scattering angle of the different energy components in the incident beam [ Fig. 5(d)]. The features in the simulated image, however, are not as wide as in the measured data [ Fig. 5(a)]. This is likely explained by the incoherent point spread of the detector.
The vertical divergence of the condensed line beam is large compared with the intrinsic width of the dynamical rocking curve of the diamond sample. The crystal therefore acts as an analyzer when rocked in the condensed line beam (see Fig. 6). The width of the rocking curve is largely determined by the divergence of the incident beam, but the finite-energy bandwidth blurs the sharp cut-off caused by the condenser lens' physical aperture. The asymmetry of the measured rocking curve [ Fig. 6(c)] reveals a misalignment of the condenser lens.
High-frequency defects in the condenser lens are visible in the spatially resolved rocking curve of perfect parts of the crystal [see Fig. 6(b)] as vertical stripes. The stripes are blurred along the rocking-angle direction due to the finite bandwidth of the incident X-rays -an effect which is confirmed by the simulations shown in Fig. 6(a).
When the crystal is rocked, a different part of the spectrum of the incident beam will be in the Bragg condition and therefore the sample will scatter in slightly different directions as a function of the rocking angle [sketched in Fig. 6(d)]. Due to the finite propagation distance from the sample to the detector, this change in angle translates to a change in position of the measured image. This mix of position and angular information, which is avoided by the use of an objective lens, is unavoidable in X-ray topography methods due to the finite propagation length, but it is exaggerated in this study due to the relatively large propagation distance and large vertical divergence compared with more usual topography experiments.

DFXM images
To simulate the DFXM images, we use a model of an ideal CRL. The physical aperture of the CRL is defined by a 0.1 Â 0.1 mm square absorbing slit placed at the exit of the lens. The only fitted parameters for the whole simulation are the relative positions of the sample, lens and detector as well as the noise level of the detector. The intensity of the simulated images is scaled to match the intensity of the measured images. The same scaling parameter is used for all images. The experimental images are overexposed (saturated) at the direct image of the stacking fault; therefore we here choose a color map that clips the highest intensities in the simulated images.
As can be seen in Fig. 7, the DFXM simulation qualitatively recreates the features of the experimental images. However, there are a number of deviations:  (i) We underestimate the magnification of the imaging setup by about 5%, which leads to an incorrect scaling of the images. This is most likely due to a small deviation of the apex radius of curvature from the nominal value of 50 mm in the individual lenses of the CRL used as objective lens.
(ii) The simulated images contain a smaller number of Borrmann fringes (the horizontal stripe features) than the measured data. We attribute this to the known high sensitivity of the spacing of Borrmann fringes to small macroscopic strain gradients (Rodriguez-Fernandez et al., 2021). Alterntively, a slight miscalibration of the photon energy or sample thickness could also lead to this difference.
(iii) The simulated images have a regular pattern of vertical streaks close to the right-hand side of the images. These are due to Fresnel diffraction from the hard edge of the square aperture in the objective lens. This is likely an artifact of the assumption of perfect transverse coherence in the horizontal direction or of the perfectly sharp edges of the aperture that are somewhat jagged in practice.
(iv) The measured images contain noise with the appearance of vertical streaks and speckle-like features close to the brightest features. This can be explained either by the aberrations in the condenser lens or in the objective lens, as will be discussed later.
In Figs. 7(c), 7(d) we simulate an image where the objective lens is displaced from the center of the scattered beam such that rays that are specularly reflected fall outside the aperture of the objective lens in the bottom part of the displayed ROI (region of interest). In that region, only diffusely scattered X-rays will contribute to the image. This results in the disappearance of the dynamical features, while the direct image of the stacking fault can still be seen. In visible light microscopy, this is referred to as 'dark-field contrast'.
In the geometric optics treatment, weak beam contrast is explained by the presence of small regions where the lattice is strained and rotated away from the average lattice. In these regions, rays are scattered if they satisfy the exact Bragg condition for the deformed lattice. A stacking fault is in principle a perfectly sharp defect with no spatial extent so no such region exists. The appearance of weak beam contrast therefore illustrates the inability of the geometric model to handle diffraction effects that are important when describing scattering from small structures with a characteristic size on the order of the wavelength or smaller.    Since the stacking fault is thought to be a perfectly sharp defect, the width of the image of this defect can serve as a rough estimate of the resolution of the instrument. The stacking fault is a 2D feature, and therefore the width of the image is not only determined by the resolution of the imaging optics, but also by the projection of the part of the stacking fault illuminated by the sheet beam along the scattered beam direction.
In Fig. 8 (inset) we compare the width of this feature in the experimental images with that in the simulated images. We see that the polychromaticity does not contribute significantly to the width of the feature in the simulations. Previous studies of the chromatic aberrations in CRLs, using the same computational approach as we apply here, also find that the chromatic aberration only adds a small part to the point spread of the imaging optics (Pedersen et al., 2018).
The experimental image has an FWHM of 1.3 mm, about 0.5 mm wider in the demagnified sample plane coordinates than that predicted by the simulations. We believe that the resolution of the experiment is degraded by aberrations in the CRL lenses.

Aberrated lenses
So far, we have ignored the effect of the aberrations in the lenses. In transmission images (Lyatun et al., 2020) and Bragg images taken without the condenser lens, short-wavelength aberrations are known to cause strong speckle-like noise in the final images. The apparent absence of this noise in DFXM images is surprising at first. However, as previously observed, this noise is averaged out when the imaged field is only partially coherent (Falch et al., 2019;Carlsen et al., 2022). Normally we think of the dynamically scattered X-rays as highly coherent as the Bragg scattering effectively collimates the incident radiation, but this argument does not consider the polychromaticity of the incident radiation.
While CRLs have been shown to be nearly achromatic over the bandwidth of the monochromated beam (Pedersen et al., 2018), Bragg scattering is not: a higher-energy component of the incident beam scatters at a smaller angle and vice versa, as sketched in Fig. 5(d). Since the incident beam has a large vertical divergence (compared with both the energy bandwidth and the Darwin width) set by the aperture of the condenser lens, the integration over energies corresponds to integrating over a small spread of angles of the scattered beam. This integration averages out the high-frequency parts of the aberrations in the vertical direction, leaving features elongated in the vertical direction. A relative energy spread of 1.0 Â 10 À4 corresponds to an angular difference of 1.0 Â 10 À4 rad Â tan which gives 8 mm at the lens planecomparable with the average grain size (15 mm) of the O30Hgrade beryllium (Lyatun et al., 2020)

used in our CRL.
This effect is demonstrated in Figs. 9(a), 9(b), where an aberrated lens is constructed by multiplying the wavefront by an aberration function at the position of the first lens and the last lens in the CRL. The aberration functions used here are pure phase objects and are made by randomly placing a number of circles of random size. The amplitude and number of circles are chosen to make the simulated and measured images similar. The partially averaged speckle noise has the appearance of vertical stripes and is qualitatively similar to that observed in the real images [ Fig. 7(a)]. In a typical experiment, we do not acquire sufficient data to uniquely determine the aberrations, but an effective aberration function can be recovered using Fourier ptychography (Carlsen et al., 2022). The qualitative similarity between simulated and measured images confirms that the vertical stripe artifacts seen in DFXM images of highly perfect crystals can be    explained by high-frequency errors in the objective lens, which we know to be present.
In Figs. 9(c), 9(d) we investigate the effect of similar aberrations, to those used in the objective lens, in the condenser lens. Once again, averaging over the energy bandwidth significantly reduces the strength of the noise and results in vertical stripes. In this case the stripes are unbroken and can be followed from the top to the bottom of the image, in contrast to the noise observed in the real images and with an aberrated objective lens.

Conclusion
DFXM is based on well known physics and we can predict the images it will produce -if we know the structure of the sample. It is possible to simulate the full FOV of the prototypical DFXM instrument at ID06-HXM at the ESRF.
Comparison of our simulations and experimental findings from a near-perfect single-crystal diamond suggests that most deviations between our simulations and the observed images can be explained by non-ideal behavior of the lenses. This suggests that the performance of DFXM instruments is critically limited by the quality of the objective lens.
In general, we do not have a sufficiently accurate model of the sample structure to do full simulations of the DFXM experiment, and the data collected in a typical experiment are not enough to build a full 3D model of the sample at sufficiently high resolution. Nevertheless, simulations like the ones shown here should prove useful for evaluating possible upgrades of the instrumentation and to qualitatively study the type of contrast observed from different types of defects in near-perfect single crystals, such as isolated dislocations, twin boundaries and point defects.
More deformed crystals are difficult to simulate, both because models of the displacement field in such crystals are not easily obtainable and because the large strain would require impractical, small step sizes to ensure proper sampling of the scattering function. In these samples, dynamical diffraction effects are not thought to be important and the speckle-like noise due to lens aberrations should also be less strong, as the scattering is more diffuse. So a wavefront-based simulation approach is less appropriate in this type of sample. It may however be interesting to investigate the transition from the dynamical patterns from near-perfect crystallites to kinematical scattering from deformed crystals using our new approach.